;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
;
;    This file is part of ICTP RegCM.
;    
;    Use of this source code is governed by an MIT-style license that can
;    be found in the LICENSE file or at
;
;         https://opensource.org/licenses/MIT.
;
;    ICTP RegCM is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
;
;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"    
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
;*************************
; read in data
;*************************
  a = addfile("CAS_SRF.1990010100.nc","r")
  prec = a->pr

  xlat = a->xlat
  xlon = a->xlon
  dims = dimsizes(xlat)
  nlat = dims(1)                            
  nlon = dims(0)

  pdims = dimsizes(prec)
  ntimes = pdims(0)

  times = a->time;

;***********************
; plot
;***********************
  wks = gsn_open_wks("ps","RegCM_Precipitation")
  gsn_define_colormap(wks,"gui_default")
  
  res            = True                    ; plot mods desired
; !!!!! any plot of data that is on a native grid, must use the "corners"
; method of zooming in on map.
  
  res@mpLimitMode       = "Corners"            ; choose range of map
  res@mpLeftCornerLatF  = xlat(0,0)
  res@mpLeftCornerLonF  = xlon(0,0)
  res@mpRightCornerLatF = xlat(nlon-1,nlat-1)
  res@mpRightCornerLonF = xlon(nlon-1,nlat-1)

; The following 4 pieces of information are REQUIRED to properly display
; data on a native grid.  This data should be specified somewhere in the
; model itself.
  
  prj = a@projection
  clon = a@longitude_of_projection_origin

  if (prj .eq. "LAMCON") then
    trlats = a@standard_parallel
    res@mpProjection        = "LambertConformal"
    res@mpLambertParallel1F = trlats(0)
    res@mpLambertParallel2F = trlats(1)
    res@mpLambertMeridianF  = clon
  end if
  if (prj .eq. "NORMER") then
    res@mpProjection        = "Mercator"
  end if
  if (prj .eq. "POLSTR") then
    clat = a@latitude_of_projection_origin
    res@mpProjection        = "Stereographic"
    res@mpRelativeCenterLon = True
    res@mpCenterLonF = clon
    res@mpRelativeCenterLat = True
    res@mpCenterLatF = clat
  end if
  if (prj .eq. "ROTMER") then
    clat = a@latitude_of_projection_origin
    res@mpProjection        = "Mercator"
    res@mpCenterLonF = clon
    res@mpCenterLatF = clat
  end if
  
; usually, when data is placed onto a map, it is TRANSFORMED to the specified
; projection. Since this model is already on a native lambert conformal grid,
; we want to turn OFF the tranformation.
  
  res@tfDoNDCOverlay         = True      ; do not transform
  res@cnFillOn               = True      ; color plot desired
  res@cnLinesOn              = False	   ; no contour lines
  
  res@mpGeophysicalLineColor = "red"     ; color of continental outlines
  res@mpPerimOn              = True      ; draw box around map
  res@mpGridLineDashPattern  = 2         ; lat/lon lines as dashed
  
  res@gsnAddCyclic           = False     ; regional data don't add

  res@pmTickMarkDisplayMode  = "Always"  ;
  res@lbAutoManage           = True      ;

  res@tiMainFuncCode = "~" ;
  i = 1
  do while (i .lt. ntimes)
    res@tiMainString = a@title + "~C~" + ut_string(times(i), "");
    plot = gsn_csm_contour_map(wks,prec(i,:,:),res)
    i = i + 1
  end do

end
